function output = run_income_process(data,col_long,fp,output_dir)

NLI_data_mom = NLI_moment_funct(data,col_long);
output.NLI_data_mom = NLI_data_mom;

do_boot_NLI = 1; 
% Set this dummy (do_boot_NLI) equal to 1 if you want to calculate and save the weighting matrix for NLI process
% Or set it equal to 0 if you already ran this code before and want to save time
% Note: this dummy must be set equal to 1 if running the code for the first time!

% Alternative: bootstrap a weighting matrix for the NLI process estimation
%(block bootstrap over fp.N number of children)
if do_boot_NLI == 1
    display('Calculating bootstrapped data moment weight matrix for NLI process - ONLY NEED TO DO THIS ONCE');
    tic
    %pre-define
    NLI_boot_moments = zeros(length(NLI_data_mom),fp.B);

    %random number draws for bootstrap
    %this creates fp.N x fp.B matrix of random data integers between 1 and fp.N, with replacement
    rand('state',1568);
    draws = randi(fp.N,fp.N,fp.B);

    for b = 1:fp.B
        %randomly re-draw with replacement from original data
        fp.rand_order = draws(:,b); % fp.N by 1 vector of integers (with replacement)
        boot_data = [];
        for q = 1:fp.N
            %block for mother with id == rand_order(q)
            boot_data = [boot_data; data( (data(:,col_long.id) == fp.rand_order(q)),:)];
        end
        %calculate moments for data draw
        temp = NLI_moment_funct(boot_data,col_long); % using bootstrapped data now! See beginning of moment_funct.m
        NLI_boot_moments(:,b) = temp;
    end

    %which elements of moment vectors have any variation?
    %--some moments defined have NaN since out of age range, etc.
    temp = std(NLI_boot_moments,[],2);
    NLI_moments_with_variation = (temp > 0); % save the indices of the used moments
    NLI_boot_moments = NLI_boot_moments(NLI_moments_with_variation,:);
    [NLI_numb_moments,fp.B] = size(NLI_boot_moments);

    %create bootstrapped covariance matrix, off-diagonal elements set to zero
    cov = zeros(NLI_numb_moments,1);
    for b = 1:fp.B
        diff = NLI_boot_moments(:,b) - mean(NLI_boot_moments,2);
        cov = cov + (1/fp.B).*diff.^2;
    end
    NLI_wt_matrix = inv( diag(cov) );

    save([output_dir 'temp_NLI_weight_matrix'],'NLI_wt_matrix');
    save([output_dir 'temp_NLI_variation'],'NLI_moments_with_variation');

    toc
end

load([output_dir 'temp_NLI_weight_matrix'],'NLI_wt_matrix');
fp.NLI_wt_matrix = NLI_wt_matrix;

load([output_dir 'temp_NLI_variation'],'NLI_moments_with_variation');
fp.NLI_moments_with_variation = NLI_moments_with_variation;

% Now estimate the NLI parameters
estimate_NLI_process = 0;
% Set this dummy equal to 1 if you want to (re-)estimate the non-labor income process parameters
% Set it equal to 0 if you already ran this code before and want to save time
% Note: this dummy must be set equal to 1 if running the code for the first time!

if estimate_NLI_process == 1

    disp('Estimating Non-labor income process using fminsearch... - ONLY NEED TO DO THIS ONCE')
    tic
    % Load initial (or existing) estimates
    load([output_dir 'temp_estimates_NLI_process'],'NLI_par');
    param_start_NLI = NLI_par;
    % Use fminsearch to estimate NLI parameter using simulated MoM:
    [param_est_NLI,obj_eval_est,exitflag] = fminsearch_matt(@(param) NLI_obj_funct(param,fp,data,col_long,NLI_data_mom,estimate_NLI_process), param_start_NLI,fp.options_est_NLI);
    toc

    NLI_par = param_est_NLI; % these are untransformed!!
    save([output_dir 'temp_estimates_NLI_process'],'NLI_par');
end

% Note: if estimate_NLI_process = 0, we just simulate NLI data at the current set of NLI parameters
% Load most recent set of estimates
load([output_dir 'temp_estimates_NLI_process'],'NLI_par');

disp('Simulating NLI at current parameters...')
estimate_NLI_process = 0; % KEEP THIS FIXED - needed for simulation
tic
output = NLI_obj_funct(NLI_par,fp,data,col_long,NLI_data_mom,estimate_NLI_process);
toc
end % function